# -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 09:30:27 2021

@author: QCH
"""
import os
cur_path = os.getcwd()
os.chdir( cur_path )

import pandas as pd
import numpy as np
import math
import time
t_start = time.time()

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from scipy.interpolate import griddata


'set grid'
x_size = (110) *10**-3 #尺寸，单位m
y_size = (110) *10**-3
z_size = (50)  *10**-3

x_coord_start = -0.5 * x_size
y_coord_start = -0.5 * y_size
z_coord_start = 0

grid_dist = 1 *10**-3 #每个正方体cell的边长，单位m


x_num = math.ceil(x_size / grid_dist)
y_num = math.ceil(y_size / grid_dist)
z_num = math.ceil(z_size / grid_dist)
# grid_points = np.zeros((x_num*y_num*z_num,3),dtype = float)


r2 = 0.5 * x_size #半径
r1 = 0.5 * x_size * 0.98 #半径

'initialize the point in geometry flow'
shell_points = np.zeros((x_num*y_num*z_num,3),dtype = float)
shell_position = np.zeros((x_num*y_num*z_num,3),dtype = int)
grid_num = np.shape(shell_position)[0]
point_index = -1
for m in range(x_num):
    for n in range(y_num):
        for k in range(z_num):
            radius = ((x_coord_start + m*grid_dist)**2 + (y_coord_start + n*grid_dist)**2)**0.5
            x_coord = x_coord_start + m*grid_dist
            z_coord = z_coord_start + k*grid_dist
            check1 = (  radius<r2 and radius>r1 and x_coord < 0             )
            check2 = (  radius<r2 and z_coord < (z_coord_start+0.01*z_size) )

            if(check1 or check2):
                point_index = point_index + 1
                #存储xyz坐标
                shell_points[point_index] = [m*grid_dist + x_coord_start,    \
                                            n*grid_dist + y_coord_start,    \
                                            k*grid_dist + z_coord_start ]
                #存储mnk值    
                shell_position[point_index]=[m,n,k]
                
shell_points = shell_points[0:point_index]
shell_position = shell_position[0:point_index]

colname = ['x','y','z']
shell_dataframe = pd.DataFrame(data = shell_points,columns = colname)
shell_dataframe.to_csv('data\\shell_points_half.csv',index=None)


